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ABSTRACT 

We present the results from axisymmetric time-dependent hydrodynamical calculations of gas flows 
which are under the influence of gravity of a black hole in quasars. We assume that the flows are non- 
rotating and exposed to quasar radiation. We take into account X-ray heating and the radiation force due 
to electron scattering and spectral lines. To compute the radiation field, we consider an optically thick, 
geometrically thin, standard accretion disk as a source of UV photons and a spherical central object as 
a source of X-rays. The gas temperature and ionization state in the flow are calculated self-consistently 
from the photoionization and heating rate of the central object. 

We find that for a 10^ M0 black hole with an accretion luminosity of 0.6 of the Eddington luminosity 
the flow settles into a steady state and has two components: (1) an equatorial inflow and (2) a bipolar 
inflow/outflow with the outflow leaving the system along the disk rotational axis. The inflow is a 
realization of a Bondi-like accretion flow. The second component is an example of a non-radial accretion 
flow which becomes an outflow once it is pushed close to the rotational axis where thermal expansion 
and the radiation pressure accelerate it outward. Our main result is that the existence of the above 
two flow components is robust to the outer boundary conditions and the geometry and spectral energy 
distribution of the radiation field. However, the flow properties are not robust. In particular, the 
outflow power and coUimation is higher for the radiation dominated by the UV/disk emission than for 
the radiation dominated by the X-ray/central engine emission. Our most intriguing result is that a 
very narrow outflow driven by radiation pressure on lines can carry more energy and mass than a broad 
outflow driven by thermal expansion. 

Subject headings: accretion, accretion disks - methods: numerical - HD 



1. INTRODUCTION 

Active Galactic Nuclei (AGN) in general, but quasars 
in particular are very powerful sources of radiation. The 
AGN luminosity, L is typically high compared to its Ed- 
dington hmit, LEdd (0.001 ^Edd ^ L ^1 ^Edd)- The spec- 
tral energy distribution (SED) of AGN is very broad. It 
spans the wavelength range from radio to hard X-rays and 
even TeV. Most of the AGN luminosity is in the optical- 
UV-IR bands but some significant fraction is in the X-ray 
band. 

These AGN radiation properties and the AGN central 
location in their host galaxies imply that they play a very 
important role in determining the ionization structure and 
dynamics of matter not only in their vicinity but also 
on larger, galactic and even intergalactic scales (Ciotti 
& Ostriker, 1997, 2001; King 2003; Murray, Quataert, & 
Thompson 2005; Sazonov et al. 2005; Springel, Di Matteo 
& Hernquist 2005; Hopkins et al. 2005; Wang, Chen, & 
Hu and references therein). There are many indications 
that support this suggestion, for example the presence of 
broad emission lines (BELs) in AGN spectra. 

BELs are one of the defining spectral features of AGN. 
They are observed in optical (OPT) and ultraviolet (UV) 
spectra and have line wings extending to velocities up to 
10^ km s^^. It is fairly well established that the primary 
physical mechanism for production of BELs is photoion- 
ization by the compact continuum source of AGN. De- 
tailed photoionization calculations yield relatively tight 
constraints on the physical conditions of the emitting gas 
- such as the temperature, density, ionization level, and 



chemical abundances (e.g., Ferland et al. 1998; Hamman 
& Ferland 1999; Krolik 1999 and references therein). The 
width of BELs indicates that the photoionized gas is su- 
personic. The shape and position of the line profiles can 
be explained by lines forming either in a region without 
preferred velocity direction and with a nearly spherical 
distribution or at the base of a wind from an accretion 
disk (Murray et al. 1995, hereafter MCGV; Bottorff et. al 

1997) . In the latter case, the BELs can form where the 
outflow expansion velocity is low compared to the rota- 
tional velocity (MCGV). 

Not all gas in AGN shares the properties of the one 
responsible for BELs. Some quasars show broad absorp- 
tion lines (BALs) which are the most dramatic evidence 
for well-organized outflows in AGN. BALs are almost al- 
ways blueshifted relative to the emission-line rest frame, 
indicating the presence of outflows from the active nu- 
cleus, with velocities as large as 0.2 c (e.g., Turnshek 

1998) . BALs are observed not only in the UV but also 
in the X-rays. For example, Chartas, Brandt & Gallagher 
(2003) discovered a very broad absorption line in the X-ray 
spectrum of PG 1115-1-80. Other evidence for AGN out- 
flows include narrow absorption lines (NALs). UV spectra 
of some quasars show NALs which are blueshifted by as 
much as ~ 50000 km s~^ (Hamann et al. 1997). NALs are 
found much more commonly in the UV spectra of Seyfert 
galaxies than of quasars but in Seyfert galaxies the lines 
are blueshifted only by several 100 km s~^ (Crenshaw et 
al. 1999). As BALs, NALs are observed not only in the 
UV but also in the X-rays. For example, Kaastra et al. 
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(2000) observed NALs due to highly ionized species in 
a high-resohition X-ray observation of the Scyfert galaxy 
NGC 5548 obtained by Chandra. As mentioned above, 
some argue that BELs may also be associated with the 
base of a disk outflow. 

The dynamics of the gas responsible for BALs, BELs, 
and NALs can be driven by radiation, even for sub- 
Eddington sources. The driving can be due to radiation 
pressure or radiation heating, or both. The radiation force 
can overcome gravity for sub-Eddington sources when the 
gas opacity is higher than the electron scattering. The 
latter is usually used to define Lsdd- The gas opacity 
can be enhanced by scattering of photons by UV spec- 
tral lines. Radiation pressure on spectral lines (line force) 
can be significant, provided gas is moderately ionized and 
can interact with the UV continuum through very many 
UV line transitions. For highly ionized gas, line force is 
negligible because of a lower concentration of ions capa- 
ble of providing UV line opacity. In the case of highly or 
fully ionized gas, an outflow still can be produced if the 
gas heating is efficient enough for the thermal energy to 
exceed the gravitational energy. AGN with their broad 
SED, are systems where both line driving and radiation 
heating, in particular. X-ray heating can operate. 

The gas dynamics can be also affected by dust. For 
example radiation pressure due to the dust opacity can 
produce an outflow of dust and also of g the two 

can be coupled (e.g., Phinney 1989; Pier & Krolik 1992; 
Emmcring, Blandford, & Shlosman 1992; Laor & Draine 
1993; Konigl & Kartje 1994; Murray et al. 2005 and refer- 
ences therein). Matter in AGN is known to be a mixture 
of gas and dust on large scales (e.g.. Miller & Goodrich 
1985; Awaki et al. 1991; Antonucci 1984; Blanco, Ward, 
& Wright 1990; Krolik 1999 and references therein) but 
also on subparsec scales as in the Seyfert galaxy NGC 
1068 (e.g., Wittkowski et al. 2004). However, details of 
the matter composition and microphysics are very com- 
plex and poorly understood. Additionally, the presence of 
dust affects the state and dynamics of gas and vice verse. 
In particular, it is nontrivial to include the dust effects on 
the gas photoionization state. Therefore studies that con- 
sider both matter components rely on simplifying assump- 
tions regarding the dynamics and geometry (e.g., Murray, 
Quataert, & Thompson 2005). There are also studies that 
focus only on the gas component and consider dynamics 
on relatively smaller radii (e.g., Sazonov et al. 2005). We 
follow the latter and focus on the gas dynamics. 

It is commonly accepted that AGN are powered by disk 
accretion of matter onto a supermassive black hole (BH). 
A wind driven from such a disk by line force is the most 
promising hydrodynamical (HD) scenario for AGN out- 
flows, especially high luminsity quasars. In this scenario, 
a wind is launched from the disk by the local disk radia- 
tion at radii where the disk radiation is mostly in the UV 
(Shlosman, Vitello & Shaviv 1985; MCGV). Such a wind 
is continuous and has mass loss rate and velocity which 
are capable of explaining the blueshifted absorption lines 
observed in many AGN, if the ionization state is suitable 
(e.g., MCGV, Proga, Stone & Kallman 2000, PSK here- 
after; Proga & Kallman 2004). This wind scenario has the 
desirable feature that for the wind motive power, it relies 
on radiation which is an observable quantity. However, 
not all AGN outflows can be explained by line driving be- 



cause of too low himinosity or too high ionization state or 
both (e.g.. Chelouche & Netzer 2005; and Kraemer et al. 
2006). Therefore, other mechanisms such as thermal and 
magnetic driving are likely important. 

Theoretical models predict that X-ray heating can have 
profound effects on the gas dynamics in disks. Since X-rays 
tend to heat low density gas to a temperature Tc ~ lO'' K, 
with which matter in an accretion disk is expected to ei- 
ther puff up and form a static corona, or to produce a 
thermal wind, depending on whether the thermal veloc- 
ity exceeds the local escape velocity, v^sc (e.g. Begelman, 
McKee and Shields, 1982; Ostriker, McKee, & Klein 1991; 
Woods et al. 1996; Proga & Kallman 2002). This paper, 
and this section in particular, are focused on non-magnetic 
processes. Magnetic driving is briefly discussed in §4 (see 
also Proga 2007). 

The insights gained from these studies support the radi- 
ation driven disk wind model for outflows in quasars and 
likely in other AGN. These successes motivate further ex- 
ploration of radiation effects on gas dynamic's in and out- 
side AGN using similar methods. The particular issues 
needed attention include: (1) Can AGN radiation drive 
an outflow from anywhere else than an accretion disk (e.g., 
from a large scale torus believed to surround an AGN or 
from an inflowing gas accreting directly onto a BH or in- 
directly through an accretion disk)? (2) If so, what is the 
power of such an outflow and can AGN radiation control 
the rate at which matter is supplied to the AGN accretion 
disk. More broadly, is the AGN radiation an important 
element of the so-called AGN feedback? (3) What is the 
origin and the dynamics of matter responsible for NALs 
and NELs? 

We report here on results from our first phase of study- 
ing gas dynamics in AGN on sub- and parsec-scales. We 
have simplified several aspects of this complex problem. 
For example, we have not taken into account rotation of 
the fiow and the magnetic fields. Additionally, we have not 
taken into ac:count dust as we assume that all the dust has 
been destroyed. This assumption is valid interior to the 
dust sublimation radius, where the central object radia- 
tion is responsible for the dust destruction. At large radii, 
still inside our computational domain, this assumption is 
less valid but the dust could be destroyed there by other 
processes such as dust sputtering. Our goal is to set-up 
simulations as simple as possible to study the effects of 
the X-ray heating (important in the so-called preheated 
accretion, sec Ostriker et al. 1976) and radiation pressure 
on gas which is initially captured by a BH. 

We decided to study non-rotating gas because once this 
gas is initially located in the polar region at large radii, 
it will likely stay in this region regardless of the subse- 
quently radial movement with respect to the BH and its 
accretion disk, i.e, it will not settle down on the accretion 
disk. Therefore non-rotating gas can be exposed to the 
direct disk radiation for a wide range of radii and radial 
velocities. For these reasons, non-rotating gas can easily 
absorb, re-emit, and scatter the disk radiation and be re- 
sponsible for line emission and absorption in AGN. 

We use the methods developed by PSK to study gas 
dynamics on sub-parsec- and parsec-scales in AGNs. We 
consider an axisymmetric HD flow accreting onto a super- 
massive BH. The flow is non-spherical because it is irradi- 
ated by an accretion disk and spherical corona. The disk 
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radiation flux is the highest along the disk rotational axis 
and is gradually decreasing with increasing polar angle, 
9 as cos9. The corona radiation is isotropic. We take 
into account the radiation heating and cooling, radiation 
pressure due to the electron scattering and spectral lines. 
We adopt a simplified treatment of photoionization, and 
radiative cooling and heating that allow us to compute 
self-consistently the ionization state, and therefore the line 
force, in the flow. 

We concentrate on assessing what fraction of the inflow 
from large radii reaches the vicinity of a BH and what 
fraction is turned into an outflow. We also assess the im- 
pact of thermal expansion and radiation pressure on driv- 
ing the outflow and on the outflow properties such as the 
two-dimensional structure, thermal and kinetic power and 
temporal behavior. We describe our calculation in Section 
2. We present our results in Section 3. The paper ends, in 
section 4, with discussion and our conclusions. 

2. METHOD 

2.1. Hydrodynamics 

To compute the structure and evolution of a flow irradi- 
ated by an AGN, we solve the equations of hydrodynamics 

g+pV.v = 0, (1) 
P^ = -VP + pg + pF-d (2) 

where p is the mass density, P is the gas pressure, v is the 
velocity, e is the internal energy density, £ is the net cool- 
ing rate, g is the gravitational acceleration of the central 
object, F'''^'^ is the total radiation force per unit mass. We 
adopt an adiabatic crjuation of state P = (7 — l)e, and 
consider models with the adiabatic index, 7 = 5/3. 

To solve eqs. 1-3, we use the ZEUS-2D code (Stone & 
Norman 1992) extended by PSK. We perform our calcula- 
tions in spherical polar coordinates (r, 9, <p) assuming axial 
symmetry about the rotational axis of the accretion disk 
{9 = 0"). 

Our standard computational domain is defined to oc- 
cupy the angular range 0° < 9 < 90° and the radial range 
Ti = 500 r* < r < Vq = 2.5 x 10^ r*, where r* = 3rs 
is the inner radius of the disk around a Schwarzschild BH 
with a mass, Mbh and radius rg = 2GMbh/c^ ■ The r — 9 
domain is discretized into zones. For Tq = 2.5 x 10^ 
and 1 X 10® r*, our numerical resolution in the r direction 
consists of 140 and 180 zones, respectively. We fix the zone 
size ratio, drk+i/drk = 1.04 (i.e., the zone spacing is in- 
creasing outward). Gridding in this manner ensures good 
spatial resolution close to the inner boundary at n. In 
the 9 direction, our standard numerical resolution consists 
of 50 but we also performed simulations with 100 zones. 
Zone size ratio is always dBi/dOi+i = 1.0 (i.e., grid points 
are equally spaced). 

For the initial condition, we typically assume spheri- 
cal symmetry set all HD variables to constant values ev- 
erywhere in the computational domain, i.e., p{r,9) = pQ 
e{r,9) = eo, Vr{r,9) = Vrfi = 0, ve{r,9) = vefi = 0, and 



v,p{7\9) = f0,o = 0. For comparison, we performed simu- 
lations where we initially set HD variables to values which 
resulted from one-dimensional simulations without radia- 
tion pressure. 

We specify the boundary conditions in the following way. 
At the pole (i.e., 9 = 0°), we apply an axis-of-symmetry 
boundary condition while for 9 = 90°, we apply reflect- 
ing boundary conditions. For the inner and outer radial 
boundaries, we apply an outflow boundary condition (i.e., 
to extrapolate the flow beyond the boundary, we set val- 
ues of variables in the ghost zones equal to the values in 
the corresponding active zones, see Stone & Norman 1992 
for more details). To represent steady conditions at the 
outer radial boundary, during the evolution of each model 
we continue to apply the constraints that in the last zone 
in the radial direction, v${ro,9) = v^{ro,9) = 0, and the 
density and internal energy are fixed at constant values 
at all times. We allow Vr to float. We have found that 
this technique produces a solution that relaxes to a steady 
state solution for both spherically symmetric accretion and 
non spherical accretion with an outflow (see also Proga & 
Begelman 2003a). Our approach then is to mimic the sit- 
uation where there is always gas available for accretion. 

2.2. The radiation field and force 

The geometry and assumptions needed to compute the 

radiation field from the disk and central object are sim- 
ilar to those considered by PSK but here we make some 
simplifications to the earlier model. The disk is flat, Kep- 
lerian, geometrically-thin and optically-thick. The disk 
photosphere is coincident with the 9 = 90" axis. We 
specify the radiation field of the disk by assuming that 
the temperature follows the radial profile of the opti- 
cally thick accretion disk (Shakura & Sunyaev 1973), and 
therefore depends on the mass accretion rate in the disk, 
Md, the mass of the BH, Mbh and the inner edge of 
the disk, r*. In particular, the total accretion luminosity, 
L = 2r]GMBHMs,/rs, where rj is the rest mass conversion 
efficiency. 

The geometry of the central engine in AGNs is poorly 
known. As in PSK, we consider the central engine as the 
most inner part of the accretion disk plus an extended 
corona. We refer to the corona as the central object. The 
inner radius of the computational domain is large com- 
pared to the radius of the central object. Therefore the 
central object can be approximated as a point source lo- 
cated at r = 0. We express the disk luminosity, Lu and the 
central object luminosity, in units of the total accretion 
luminosity, i.e., Ld = foL and L, = /*L = (1 — /d)^. For 
simplicity, we assume that the disk emitts only in the UV, 
whereas the central object emitts only the X-rays, i.e., 
the system UV luminosity, Luv = fwL = Ld and the 
system X-ray luminosity, Lx = fxL = L* (in other words 
/uv = /d and fx = /*)• 

With the above simplifications, only the central object 
radiation is responsible for ionizing the flow to a very high 
state. The central radiation contributes to the radiation 
force due to electron scattering but does not contribute to 
line driving. On the other hand, the disk radiation con- 
tribute to the radiation force due to both electrons and 
lines. Note that we assume that the luminosity in the re- 
maining bands, mainly optical and infrared, is the part of 
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the luminosity that docs not change the dynamics of the 
wind and set it to zc;ro. 

To evaluate, line force, we generally follow PSK who 
used a modified Castor, Abbott & Klein method (Castor, 
Abbott, & Klein 1975; see also Proga, Stone, Drew 1998, 
1999). The line force at a point defined by the position 
vector r is 
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(r) 
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where / is the frequency-integrated continuum intensity 
in the direction defined by the unit vector ri, and f2 is the 
solid angle subtended by the disk and central object at the 
point. The term in brackets is the electron-scattering ra- 
diation force, (Te is the mass-scattering coefficient for free 
electrons, and Af (t) is the force multiplier (Castor ct al. 
1975) - the numerical factor which parametrizes by how 
much spectral lines increase the scattering coefficient. In 
the Sobolev approximation, M{t) is a function of the op- 
tical depth parameter 
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where t;th is the thermal velocity, and ^ is the velocity 
gradient along the line of sight, n. 

Evaluation of the radiation force due to an non- 
isothermal extended source such as an accretion disk is rel- 
atively expensive (see PSK, Proga et al. 1999). However, 
we consider here the flow relatively far from both the cen- 
tral object and the most luminous part of the disk. There- 
fore, we simplify the force evaluation as follows. First, we 
consider only the radial component of the force, FJl''^'' and 
set the other components to zero. Second, we compute 
the velocity gradient by taking into account only the dom- 
inant term: gradient of the radial velocity along the radial 
direction, i.e. 

' dr 

With these simplifications and an assumption that the 
flow does not change the geometry of the radiation field, 
we can evaluate the radiation force in two steps: (1) cal- 
culation of the velocity gradient in the fir direction and 
then the optical depth parameter t; (2) calculation of the 
parameters of the force multiplier using a current value 
of the photoionization parameter, ^ adopting results of 
Stevens & Kallman (1990). Then we calculate the radia- 
tion force exerted by radiation along the radial direction. 
Our approach to evaluate the radiation forces is simplified 
compared to the approach used in PSK because we do not 
integrate the radiation force over the solid angle subtended 
by the radiant surface but instead use an analytic approx- 
imation (compare eqs. 4 and 8), and we do not correct the 
radiation force for the optical depth effects. 

Below we describe in more detail, our calculations of the 
force multiplier and the radiation force for various condi- 
tions in the fiow. For r ^ r*, in the optically thin case, 
the radial radiation flux from a disk and central object can 
be written as 



TTy{r,e) = 2 cos 6* 
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and 



(6) 
(7) 



respectively. Eq. 6 approximates the numerical evaluation 
of the radiation flux integral over the solid angle subtended 
by the disk surface (see eqs B2 and B3 in Proga, Stone & 
Drew 1998) while eq. 7 is the exact solution for the radi- 
ation flux for a point source. With eqs. 6 and 7, one can 
show that eq. 4 reduces to 
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(r,^) = ^^[/* + 2cose/D(l + M(t))] (8) 

To compute the force multiplier, we adopt Castor et al. 's 
(1975) analytical formula modifled by Owocki, Castor & 
Rybicki (1988 see also Proga et al. 1998) 



M{t) = kt~ 
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where k is proportional to the total number of lines, 

a is the ratio of optically-thick to optically-thin lines, 
''max — t r/max and rymax is a parameter determining the 
maximum value. Mmax achieved for the force multiplier. 
Equation 9 shows the following limiting behavior: 



lim M{t) 
lim M{t) 



= M„ 



(10) 

(11) 



where M^ax = k{l - a)ri^^. 

To compute, the parameters of M{t) taking into account 
the effects of ionization, we follow PSK. In particular, we 
compute the parameters of the line force using a value of 
the photoionization parameter, ^ and the analytical for- 
mulae for k and Tjmax from Stevens & Kallman (1990): 



k = 0.03 -I- 0.385 exp(-1.4 ^ ), 



(12) 



and 



logio Vn 



6.9 exp(0.16 C° '^) 
9.1 exp(-7.96 X 10-3^) for log 



for log;^Q ^ < 0.5 
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C > 
(13) 

The parameter a. = 0.6 and does not change with ^. 

Our procedure to calculate M{t) is computationally effi- 
cient and gives approximate estimates for the number and 
opacity distribution of spectral lines for a given ^ without 
detail information about the ionization state. However, 
it is just an estimate because Stevens & Kallman (1990) 
results were computed for the case where a 25 000 K stel- 
lar atmosphere was irradiated by a 10 keV bremsstrahlung 
ionizing spectrum whereas we consider gas not of a stellar 
atmosphere and its temperature, without irradiation, de- 
pends of the outer boimdary conditions. Additionally, the 
AGN spectrum is different from a bremsstrahlung spec- 
trum. As we describe below, the parameters k depends on 
the gas temperature (i.e., eq. 17). Additionally this and 
other parameters of M(t) depend also on the SED of the 
ionizing and UV radiation. In particular, log^g ^ = —0.5 
and the UV emission caused by the low energy extension 
of the 10 keV bremsstrahlung, Mmax ~ 300 whereas for 
an ACN-like UV spectrum, i.e., a power law with the en- 
ergy index equal to 1, M^ax ~ 8000 (T. Kallman private 
communication). For comparison, eqs. 12 and 13 yield 
Afmax ~ 3400 also for log;^o^ = —0.5. A detail investiga- 
tion of these and other dependencies is beyond the scope 
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of this study but it is important to carry out such a study 
in the future. 

To calculate the photoionization parameter needed in 
eqs. 12 and 13 and in evaluation of the net radiative cool- 
ing, we use the formula: 



47rJFx 



(14) 



where Tx is the local X-ray flux, n is the number density of 
the gas (= p/{mpy) , where nip is the proton mass, and /x 
is the mean molecular weight). We consider models with 
11 = 1. The local X-ray flux is corrected for the optical 
depth effects: 

= J^* exp(-rx), (15) 

where tx is the X-ray optical depth. We estimate rx be- 
tween the central source and a point in a flow from; 



rx 



Jo 



KxP dr, 



(16) 



where kx is the absorption coefficient, and r is the distance 
from the central source. 

The flow physical conditions and consequently the line 
force depend not only on ^. In particular, the gas tem- 
perature and ionization state can be affected by adiabatic 
heating or cooling. For example, the line force can be 
neglected even in the region of low ^ if this region is hot 
due to adiabatic compression. Therefore, we first compute 
the line force parameters based on ^ and then correct the 
parameters for the gas temperature effects. 

Specifically, we compute the parameter k using the ex- 
pression 

C -0.383 for logio T <^ 

\og-^ok={ -0.630 logio r + 2.138 for4 < logip 7" < 4.75 
I -3.870 logio T + 17.528 for logio T > 4-75 

We adopt this value of k if it is smaller than the one ob- 
tained by using expression 12. Expression 17 is based on 
detailed photoionization calculations performed using the 
XSTAR code (T. Kallman private communication). 

These expressions for the parameters of the force mul- 
tiplier predict that Mmax increases gradually from ~ 2000 
to 5000 as ^ increases from to ~ 3 and then drops to 
~ 0.1 at ^ = 1000. The line force becomes negligible for 
^ ^ 100 because then Mmax ^ 1- The line force becomes 
also negligible if T > 10^ K for any ^. 

2.3. Radiation heating and cooling 

To calculate the gas temperature we also follow PSK, 
who assumed that the gas is optically thin to its own cool- 
ing radiation. The net cooling rate depends on the density, 
p, the temperature, T, the ionization parameter ^, and the 
characteristic temperature of the X-ray radiation Tx- We 
refer a reader to PSK and Proga & Kallman (2004) for 
more details. 

We calculate the evolution of the internal energy and 
therefore we specify the initial and outer boundary condi- 
tions for e. We do this in the following way: we specify 
the gas temperature at the outer boundary (in the last 
grid zone in the radial direction), Tq, by setting it to the 



Compton temperature, Tq 
the internal energy from 



eo = 



0.25 Tx- Then we calculate 
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3. RESULTS 

We assume the mass of the non-rotating BH, 
.^BH = 10^ Mq and the disk inner radius, 
r^, = 3 rs = 8.8 x 10^^ cm throughout this pa- 
per. We consider the rest mass conversion efficiency 
rj = 0.0833. We set to be 10^^ g s'^ (=1.6 Mq yr'^). 
These system parameters yield the accretion luminosity, 
L = 7.5 X 10*^ erg s-^(= 2 x 10^^ Lq). This luminosity 
corresponds to the Eddington number, F = 0.6. We define 
the Eddington number as L in units of the Eddington lu- 
minosity for the Schwarzschild BH, iEdd = 4:WcGMBH/cre 
[i.e., F = L/LEdd = (CTeMa)/(87rcrs)]. To determine the 
radiation field, we specify also /uv and fx- We note that 
we fix Ma during our simulations, i.e. Ma does not depend 
on the mass flux through the inner boimdary, A'/in(n) (see 
eq. 20 for the formal definition of Min(ri) and §4 for dis- 
cussion of this aspect of the calculations) . 

The parameters /uv and fx (and also /d and /*) spec- 
ifying the radiation field are our free parameters which we 
focus on exploring the most in this paper. We consider 
three cases: case A with /uv = 0.5 and fx = 0.5, case B 
with /uv = 0.8 and fx = 0.2, and case C with /uv = 0.95 
and /x = 0.05 (see Table 1 for summary of our runs). The 
spectral energy distribution of the ionizing radiation is not 
well known, our choice of values for /uv and fx is guided 
by the observational results from Zheng et al. (1997) and 
Laor et al. (1997). 

To calculate the gas temperature, we assume the tem- 
perature of the X-ray radiation, Tx = 8 x lO'' K (e.g., 
Sazonov, Ostriker & Sunyaev 2005 and references therein) 
and the line cooling parameter 6=1 (see Blondin 1994 
and PSK). The force multiplier depends only formally on 
the thermal speed, vth- Therefore to compute t, we set to 
20 km s^^, i.e., the thermal speed of a hydrogen atom at 
the temperature of 25000 K for which the parameters of 
the force multiplier were computed (Stevens & Kallman 
1990). Finally, the attenuation of the X-rays is calculated 
using Kx = 0.4 g^^ cm^ for all ^. The resulting optical 
depth corresponds to the Thomson optical depth. 

For our parameters, both the so-called Compton radius, 
Rc = GMBHpnip/kTc = 8x10^^ cm = 9 x 10^ r* and 
the Bondi radius, Rb = GMbh/c^ = 4.8 x 10^^ cm = 
5.5 x 10'' (Bondi 1952) are inside our computational 
domain. We assumed here that the gas temperature at 
infinity Too = Tc = 2 x 10^ K so that the sound speed at 
infinity, Coo = (7fcTc//xmp)^/^= 4x10'^ cm s~^ and Rc = 
7i?B- For the isothermal fiow, the Bondi accretion rate is 
Mb = 3.3 X lO^Sg s-i= 0.52 Moyr^^ We note that the 
dust sublimation radius, i?sub= 2.6x10^^ cm =3x10^ r* 
(e.g., Murray et al. 2005, eq. 41) is also inside the domain. 

We consider the situation where the BH gravity is re- 
duced by the (^-dependent radiation force. Therefore, the 
effective gravity is also 9 dependent. Using eqs. 6 and 7, 
the location where thermal energy equals the effect grav- 
itational energy (i.e., a modified Compton radius) can be 
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written as 

Rc = Rc{l-T{2cosefu + f*)]. (19) 

Figure 1 presents the results for cases A, B, and C with 
po = 10^^^ g cm~'^, Vrfl = fe.o = i^0.o=O. At the outer 
radial boundary, during the evolution of each model we 
continue to apply the constraints that the density and in- 
ternal energy (the latter is computed for To = To) are 
fixed at constant values at all times. The figure shows 
the instantaneous density and temperature distributions, 
and the poloidal velocity field of the models. Addition- 
ally, it shows i?c and the contours of the Mach num- 
ber , M = ^/v^+vj/cs, where Cg = jP/p is the 
sound speed. For all three cases, the flow settles quickly 
into a steady state (within ~ a few xlO^^ s which corre- 
spond to a few dynamical time scales at the outer radius, 
r = (rl/GMBHy^^ = 9 x 10" s ). The steady state con- 
sists of two flow components (1) an equatorial inflow and 
(2) a bipolar inflow/outflow with the outflow leaving the 
system along the pole. The outflow is coUimated by the 
infalling gas. Our calculations capture the subsonic and 
supersonic parts of both the inflow and outflow. Although 
the same components can be identified in the three cases 
their size, density and temperature, and the degree of out- 
flow collimation depend on the SED and the geometry of 
the radiation field. 

In run A, X-ray heating is the strongest and it acceler- 
ates an outflow to the maximum velocity of 700 km s~^. 
Although, the gas temperature decreases in the inflow 
down to some radius, it is not low enough for line force 
to be significant. The outflow collimation by the infall in- 
creases with increasing radius. This tendency is clearly 
seen in the other two cases. 

For run B, the outflow is strongly pushed toward the 
pole at r w 1.2 X 10^ r* and the flow pattern resembles 
a bottle neck where the density increases by a couple of 
orders of magnitude. The gas temperature drops there 
to the level below 10^ K. In this low temperature region 
line driving accelerates an outflow to the velocity up to 
4000 km s~^ which is significantly higher than the veloc- 
ity of the thermal outfiow in run A. 

The outflow collimation is very strong for run C where 
X-ray heating is the smallest. The gas is siphoned off 
within a very narrow channel along the pole. In this case, 
the fraction of the computational domain occupied by the 
inflow is the largest. The shape and size of the sonic sur- 
face reflects the fact that runs A and B are similar to each 
other, but they both differ significantly from run C. 

To quantify the properties of the inflow and in- 
flow/outflow components, we computed three radial mass 
rates as a function of radius: (1) the net rate 

^90° 

M„gt(r) = 47rrM pVrSm9d0, (20) 
(2) the inflow rate 

1-90° 

Min(r) = 47rr2 / pVrSm.6de for t;^ < 0, (21) 
and (3) the outflow rate 

/.90° 

Mo„t(r) = 47rr2 / pvr sin Odd for > 0. (22) 



We also computed the outfiow power carried out in the 
form of the kinetic energy 

Pk(r) = 27rrM pv^smOde for > (23) 
and in the form of the thermal energy 

/.90° 

Pth(r) = 47rr2 / eVrSmOdO for > 0. (24) 
Jqo 

These quantities are shown in Figs. 2, 3, and 4 for run A, 
B, and C, respectively. 

Fig. 2 shows that the outflow in case A is indeed ther- 
mal: Pk < ^th for all radii; the outflow starts at radii 
large compared to Rq, its properties are strong functions 
of radius, i.e., Mout, Pk, and Pth increase with radius. In 
particular, almost all the gas inside r < 2 x 10"* ends up 
leaving the domain through the inner boundary as there 
is no or a very weak outflow driven from small radii. 

The reduction of /*(= fx) compared to /d(= /uv) has 
a few signiflcant consequences on the outflow properties: 
in cases B and C the outflow rate is appreciable even at 
small radii (see figs. 3 and 4). In particular, in case C 
the net rate at large radii is almost an order of magnitude 
smaller than the inflow rate. The latter is almost canceled 
out by the outflow rate. The overall power of the outflow 
is <; 2 orders of magnitude higher in cases B and C than in 
case A. The relative weakness of thermal driving for lower 
/x cases is obvious when Pth is compared to Pk. In run B, 
Pth increases with increasing radius for r ^ 10^. But even 
at its maximum, Pth is 3 orders of magnitude lower than 
Pk. In run C, the kinetic dominance is even greater: Pk is 
higher than Pth by 4 orders of magnitude. 

The dependence of our results on the SED {fx and /uv) 
and the radiation geometry (/* and /d) can be explained 
by the fact that as fx and decrease the X-ray heat- 
ing also decreases whereas the radiation flux in the polar 
region increases compared to the flux at high 6 (see eq. 
8). Consequently, an increasing fraction of the inflow from 
large radii is pushed toward the polar region where radi- 
ation pressure exceeds gravity. This compressed inflow is 
siphoned off by radiation pressure. Perhaps the most in- 
triguing result is that the large scale flow pattern can be 
very misleading if one wants of assess the relative strength 
of inflow and flow. Visual inspection of Fig. 1 could lead 
to the conclusion that the outflow in run C is much less 
significant that the outflow in run A because in run C the 
outflow is very narrow and barely seen. However, the op- 
posite is true as shown in Table 1. The increase in the 
mass outflow rate with decreasing fx is due to increasing 
/uv and the increasing efficiency of radiative driving at 
/x. 

Table 1 shows that the inflow rate at Tq increases with 
decreasing fx- This trend is consistent with accretion the- 
ory, in particular with Bondi's estimate for the accretion 

rate, Mb oc Too : weaker X-ray heating results in an 
increase of the accretion rate because the gravitational 
sphere of influence over gas pressure is larger. However, a 
higher Mmin at large radii does not mean that the inflow 
rate at small radii must be higher too. We find that the 
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mass outflow rate also increases with decreasing /x . Addi- 
tionally, the efBciency of turning an inflow into an outflow 



is also sensitive to fx- 



3/4, 5/8, and 



Mout(ro)/Mi,(ro) 

8/9, for /x= 0.5, 0.2, and 0.05, respectively. As we men- 
tioned above, the decrease from high to low values of fx 

corresponds to a transition from a thermally to radiation 
driven outflow. We conclude then that thermal driving is 
less efficient than radiation driving. 

One of the parameters specifying our problem is the gas 
temperature at the outer boundary, Tq. For the runs A, B, 
and C, wc set the value of To equal to the Compton tem- 
perature. This value of the gas temperature is expected if 
the gas is optically thin and heating by X-rays from the 
central object. However, it is possible that the gas at large 
radii is heated also by other processes (e.g., shocks) or the 
X-ray heated is reduced by optical depth effects. There- 
fore, To can be decoupled from Tx and be lower or higher 
than TrmC- To check how are results depend on the as- 
sumed value of To, wc rerun model B with To set to 1/10, 
1/3, and 3 of Tq, runs Bl, B2, and B3, respectively. 

As shown in Table 1, Min(ro) and Mnet(''i) increase with 
increasing To. However, one would expect the opposite 

3/2 

trend because accretion theory predicts that Mb c>c Too 
(To is a proxy of Too in the simulations). This discrepancy 
is due to the fact that the above scaling law is derived for 
the Bondi accretion problem where the gas is adiabatic or 
isothermal, i.e., where the gas temperature at the critical 
point increases with Too. However, in our problem, the 
gas is not adiabatic as we allow for radiative cooling and 
heating. This has an important effect on the gas proper- 
ties (e.g., Zel'dovich & Novikov 1971) . In particular, for 
Tq > Tq, radiative processes lead to cooling as the gas 
accretes. Consequently, the sonic radius is smaller than 
that for To = Tc- Similarly, for Tq < Tq, radiative pro- 
cesses lead to heating as the gas accretes. Consequently, 
the sonic radius is larger than that for To = Tq. There- 
fore, our results for various To reflect the basic physical 
property of accretion flows: flows with a higher temper- 
ature at the sonic radius have lower accretion rates than 
flows with a lower temperature because the gravitational 
sphere of influence over gas pressure grows with decreasing 
temperature at the sonic radius. 

4. DISCUSSION AND CONCLUSIONS 

We have calculated a scries of models for non-rotating 
flows which are under the influence of super massive BH 
gravity and radiation from an accretion disk surround- 
ing the BH. Our numerical approach allows for the self- 
consistent determination of whether the flow is gravita- 
tionally captured by the; BH or driven away by thermal ex- 
pansion or radiation pressure. We find that the flow settles 
quickly into a steady state and has two components (1) an 
equatorial inflow and (2) a bipolar inflow/outflow with the 
outflow leaving the system along the pole. The first com- 
ponents is a realization of Bondi-like accretion flow. The 
second component is an example of a non-radial accretion 
flow which becomes an outflow once it is pushed close to 
the rotational axis of the disk where thermal expansion 
and radiation pressure can accelerate the flow outward. 
Our main result is that the existence of the above two flow 
components is robust yet their properties are sensitive to 
the geometry and SED of the the radiation field and the 



outer boundaries. In particular, the outflow power and the 
degree of coUimation is higher for the model with the radi- 
ation dominated by the UV/disk emission (run C) than for 
model with the radiation dominated by the X-ray/central 
engine emission (run A). This sensitivity is related to the 
fact that thermal expansion compared to radiation pres- 
sure drives a weaker, less-coUimatcd outflow. 

Our results will likely change if we allow for signiflcant 
gas rotation. In particular, the gas will tend to converge 
toward the equator due to the combination of the cen- 
trifugal and gravitational forces. This, in turn, will likely 
weaken the outflow in the polar region because less gas 
will be pushed toward the polar region. Adding rotation 
is indeed an important next step. However, within our 
framework, it requires the introduction of additional free 
parameters. In particular, we would need to prescribe the 
rotational rate at the outer boundary. I plan to report 
on our results from some simulations with rotation in a 
companion paper. 

The treatment of the inner and outer boundaries can 
affect the results of these simulations. We set up the outer 
boundary in the spirit of representing steady conditions. 
However, in real systems but also in our simulations, the 
conditions at large radii are complex. In particular, we fix 
the density and internal energy at To at all times. This 
produces self-consistent results for accretion flows but not 
for outflows because the outflow density and internal en- 
ergy at To do not have to be as those we set there through 
implementation of the outer boundary conditions. Reruns 
of our simulations on a larger computational domain and 
with the outer boundary conditions as predicted by sim- 
ulations of galaxy evolution such as those performed by 
Springel, Di Matteo & Hernquist (2005) would be very 
important. 

It is also important to explore the effects of dust on our 
results. As we mentioned in §3, the dust sublimation ra- 
dius is inside our computational domain. The presence 
of the dust at large radii can rc(hic;e the matter density 
there, because the radiation pressure on dust can acceler- 
ate the matter outward. Consequently, the mass supply 
rate and density at small radii can decrease. This in turn 
can lead to higher gas temperature at small radii. If the 
dust reduces the gas temperature at large radii, based on 
our simulation with various Tq, we expect similar changes. 
In summary, we expect that inclusion of dust will lead to 
a lower net mass accretion rate and higher temperature at 
small radii compared to the results presented here. 

Our treatment of the inner boundary conditions is one of 
the simplest possible: free outflow - no matter enters the 
computational domain through the inner boundary. How- 
ever, as we mentioned in §1, quasars have winds which 
are likely driven from relatively small radii. Therefore, 
one would need to consider matter entering our compu- 
tational domain through the inner boundary. Even if the 
small scale quasar winds do not flow very far from the 
center, their presence could affect our results. We allow 
the radiation field from the disk and central object to be 
attenuated by the gas inside the computational domain. 
However, PSK showed that disk winds in AGN can be op- 
tically thick and they can change the geometry and SED 
of the disk and central object radiation at large radii, i.e., 
interior to the inner radius of our computational domain. 

We have performed our simulations in the HD limit and 
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have not accounted for the effects of magnetic fields. The 
dynamics of accretion flows can be significantly affected 
for magnetic fields. In particular, outflows can be mag- 
netically driven from accretion disks (e.g., Blandford & 
Payne 1982; Pelletier & Pudritz 1992 Konigl 1993; Emmer- 
ing, Blandford & Shlosman 1992; Contopoulos & Lovelace 
1994; Konigl & Kartje 1994; Contopoulos 1995; de Kool 
& Begelman 1995; Bottorff et al. 1997; Bottorff, Korista 
& Shlosman 2000; Everett, Konigl, & Arav 2002; Proga 
2003; Everett 2005). In fact, most magnetohydrodynamic 
simulations of accretion flows show outflows (e.g., Uchida 
& Shibata 1985; Stone & Norman 1994; Ouyed & Pudritz 
1997; Hawley & Balbus 2002; Proga & Begelman 2003b; 
De Villier et al. 2004; McKinney & Gammie 2004, and ref- 
erences therein) . Additionally, an outflow can be produced 
even from a flow with angular momentum so low that, if 
not for the efl[ects of magnetic fields, it would accrete di- 
rectly onto a BH without forming a disk (Proga 2005). 
Thus, magnetic fields can directly and indirectly change 
our results, e.g., magnetic disk winds, similarly to radia- 
tion driven disk winds, can change the SED and they can 
interact with other flow components. We expect that flow 
interaction can lead to time-dependent evolution and for- 
mation of inhomogeneous inflows and outflows (e.g., Proga 
& Begelman 2003b) . Therefore, it would be a natural step 
to repeat simulations similar to ours but in the magne- 
tohydrodynamical limit. We have already set the stage 
for such a study because simulations performed by Proga 
(2005) are a close magnetohydrodynamical analog of those 
presented here. We note that despite two different physical 
regimes explored here and in our previous work, similarly 
simple inflow/outflows solutions were obtained. One of the 
main motivation for this work, is to gain more insights to 
the general problem of accretion and outflow production 
where many complex time-dependent flow components are 
expected. It is also desirable to connect numerical simu- 
lations with theoretical predictions such as those obtained 
by Blandford & Begelman (1999) and Henriksen & Valls- 
Gabaud (1994). 

The SED and consequently gas temperature can be dif- 
ferent from those we have in our simulations, even in an op- 
tically thin case. One of our simplifications is an assump- 
tion that the radiation temperature, Tx (or more gener- 
ally the SED) does not change with 9. This assumption 
is motivated by X-ray observations showing that quasar 
radiation heats a low-density gas, on scales comparable to 
the Bondi radius, up to an equilibrium Compton temper- 
ature of about 2 X 10'^ K (Sazonov et al. 2005; Allen et 
al. 2006 and references therein). However, the SED and 
Tx can change with 6. For example, we assume that the 
X-ray/central object radiation does not change with 9 (see 
eq. 7) whereas the UV/disk radiation decreases with in- 
creasing 9 (see eq. 6) . Consequently one could expect the 
ratio between the X-ray and UV fitix to increase with in- 
creasing 9. We plan in the near future, to test how impor- 
tant this effect is and try to reconcile the above theoretical 
expectation with observations. 

We assume that a quasar is powered by accretion, but 
we do not couple the disk accretion rate, that determines 
radiation, to the mass accretion rate through the inner 
boundary. This Ma-Min(ri) decoupling is physically mo- 
tivated because an accretion disk is build from rotating 



gas and non-rotating gas will not significantly contribute 
to the disk formation and disk radiation. In this sense we 
study radiatively inefficient fiows accreting onto an object 
with a radiatively eflicient accretion disk. Additionally, 
we note that an accretion disk, which is located interior 
to the inner radius of the computational domain, radiates 
steadily during our simulations. The fact that the fiow 
in the computational domain settles into a steady state 
within just a few dynamical time scales shows that our 
approach is reasonable. 

Outfiows in our simulations are coUimated by the in- 
flowing gas. This coUimation can lead to formation of a 
bottle-neck-like structure as we discovered in run B. There 
are two sources of concern about the collimation: is it an 
artifact of the outer boundary or of the numerical reso- 
lution. We addressed this issue by rerunning simulation 
B and C with higher resolution in the 6 direction (100 
grid points instead 50). We found that the results for the 
higher resolutions are very similar to those with the lower 
resolution counterparts. The main difference was that the 
region with T < 10^ K for r > 1.8 x 10^ r* is somewhat 
narrower (by ~ 1°) in the higher resolution run. We also 
found that the outer boundary are not responsible for the 
collimation because our test calculations with Tq four times 
larger than in the runs presented in §3 show collimation op- 
erating in a similar way and far from the outer boundary. 
Additionally, we note that in run C, collimation begins at 
r 5 X 10** (i.e., deep inside the computational domain) 
and is disconnected from the outer boundary, and for that 
matter, from the inner boundary. To study flow stability, 
in particular, stability of the outflow, it is important to 
perform fully three-dimensional simulations which will al- 
low for non-axisymmetric effects and full development of 
fluid instabilities. 

We finish by returning to the questions we raised in §1. 
Our simulations show that AGN can have a substantial 
outflow which originates from the infalling gas. Such an 
outflow can control the rate at which non-rotating mat- 
ter is supplied to the AGN central engine as its mass loss 
rate can be significantly higher than the mass inflow rate 
at small radii. For example, in run C, as little as 10% of 
the infiow at large radii reaches small radii because 90% 
of the inflow is turned into an outflow. However, even the 
power of the strongest outflow is very low compared to 
the radiation power (i.e., for run C, Pk/L = 4 x 10~^). Fi- 
nally, the gas in our simulations can be related to material 
responsible for NELs, NALs, and warm absorbers. Addi- 
tionally, the inflowing gas could be the source of the gas 
that eventually produces BALs and BELs at radii smaller 
than the inner radius of our computational domain. We 
plan to farther explore our model and relax some of our 
assumptions in order to confirm these conclusions. 
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Table 1 
Summary of results 



Run 


/d 


/. 


/uv 


/x 


To 


Mi„(ro) 


M„ot('-i) 








Pth(ro) 












2 X lO'' K 


10^^ gs-i 


10^= gs-i 


lO^'^ gs^i 


km 


lO''" erg 


10'"' org 


A 


0.5 


0.5 


0.5 


0.5 


1 


-4 


- 1 


3 


700 


2 


4 


B 


0.8 


0.2 


0.8 


0.2 


1 


-8 


- 3 


5 


4000 


100 


0.8 


Bl 


0.8 


0.2 


0.8 


0.2 


1/10 


-0.5 


- 0.09 


0.41 


1500 


0.5 


0.8 


B2 


0.8 


0.2 


0.8 


0.2 


1/3 


-2 


- 0.4 


1.6 


1700 


2 


2 


B3 


0.8 


0.2 


0.8 


0.2 


3 


-9 


- 5 


4 


400 


3 


0.8 


C 


0.95 


0.05 


0.95 


0.05 


1 


-9 


-1 


8 


6700 


300 


0.03 
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Fig. 1. — Comparison of the results for run A, B, and C (left, middle, and right column, respectively). Top row of panels: Maps of 
logarithmic density overplotted by the poloidal velocity. For clarity, the arrows are plotted with the maximum velocity set to 1000 km s^^. 
Middle row of panels: Maps of logarithmic temperature overplotted by the direction of the poloidal velocity. The solid curve in the bottom 
left corner marks the position of the Compton radius corrected for the effects of radiation pressure due to electron scattering (see eq. 19 in 
the main text). Bottom row of panels: Contours of the Mach number. The length scale is in units of the inner disk radius (i.e., r' = r/r^ 
and z' = z/rt). 
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Fig. 2. — Top panel: The mass flux rates as a function of radius for run A. The solid, dashed, and dotted line corresponds to the outflow, 
inflow, and net rates, respectively (see eqs. 20, 21, and 22 for the formal deflnitions). Note that the absolution value of the inflow and net 
rates are plotted because these quantities are negative. Bottom panel: The energy fluxes carried out by the outflow as a function of radius 
in run Al. The solid and dashed line corresponds to the kinetic and thermal energy flux, respectively (see eqs. 23 and 24 for the formal 
definitions). The length scale is in units of the inner disk radius (i.e., r' = r/rt)- 
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As Fig. 2 but for run B. 
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